/*
 * @author hzq
 * @status done
 * @see 🚀 https://www.buildenvi.com/gongju/psychrometrics
 * @date 2025-01-09 22:54:16
 * @description 湿空气状态类
 * @usage
 */

import {
  AirState,
  DryBulbTemperature,
  Moistrue,
  RelativeHumidity,
  WetBulbTemperature,
  DewPointTemperature,
  Enthalpy
} from './airState'

/**
 * ## 湿空气状态类
 */
export class Air {
  //#region 公用常量
  /**
   * 大气压力常压；单位KPa
   */
  static readonly atm = 101.325
  /** 开氏温度；单位：K */
  static readonly T0 = 273.15
  /** 干空气气体常数；单位：J/(Kg·K) */
  static readonly Ra = 287.042
  /** 水蒸气气体常数；单位：J/(Kg·K) */
  static readonly Rv = 461.524
  /** 分子量比值；molecular weight ratio */
  // static readonly MWR: 0.621945

  /** 干空气的定压比热; 单位：kJ/(kg·℃)*/
  static readonly Cpa = 1.006
  /** 水蒸气的定压比热; 单位：kJ/(kg·℃)*/
  static readonly Cpv = 1.864
  /** 水蒸气所需要的汽化潜热 */
  static readonly latentHeat = 2501
  /** 计算精度 */
  static readonly EPSILON = 1e-6
  //#endregion

  //#region 属性
  /** 大气压力属性，单位：KPa */
  pressure: number
  /** 水蒸气分压力值，单位：KPa */
  pv: number
  /** 干球温度，单位：℃ */
  t: number
  /** 含湿量，单位：Kg/Kg */
  d: number
  /** 湿球温度，单位：℃ */
  tw: number
  /** 露点温度，单位：℃ */
  td: number
  /** 相对湿度，单位：% */
  rh: number
  /** 焓值， 单位：KJ/Kg */
  h: number

  protected _pvs?: number
  //#endregion

  //#region 计算属性
  /** 湿空气密度，单位：Kg/m3 */
  get density() {
    return Air.density(this.t, this.pv, this.pressure)
  }
  /** 饱和水蒸气分压力值，单位：KPa  */
  get pvs() {
    this._pvs ??= Air.pvs(this.t)
    return this._pvs
  }
  set pvs(value: number) {
    this._pvs = value
  }
  //#endregion

  //#region 构造函数
  /**
   * ### 根据干球温度、湿球温度、大气压力构造湿空气对象
   * @param dryBulbTemperature 干球温度,单位：℃
   * @param wetBulbTemperature 湿球温度，单位：℃
   * @param pressure 大气压力，单位：KPa；默认为101.325KPa
   */
  constructor(
    dryBulbTemperature: DryBulbTemperature,
    wetBulbTemperature: WetBulbTemperature,
    pressure?: number
  )
  /**
   * ### 根据干球温度、露点温度、大气压力构造湿空气对象
   * @param dryBulbTemperature 干球温度,单位：℃
   * @param dewPointTemperature 露点温度，单位：℃
   * @param pressure 大气压力，单位：KPa；默认为101.325KPa
   */
  constructor(
    dryBulbTemperature: DryBulbTemperature,
    dewPointTemperature: DewPointTemperature,
    pressure?: number
  )
  /**
   * ### 根据干球温度、含湿量、大气压力构造湿空气对象
   * @param dryBulbTemperature 干球温度,单位：℃
   * @param moisture 含湿量，单位：Kg/Kg
   * @param pressure 大气压力，单位：KPa；默认为101.325KPa
   */
  constructor(
    dryBulbTemperature: DryBulbTemperature,
    moisture: Moistrue,
    pressure?: number
  )
  /**
   * ### 根据干球温度、相对湿度、大气压力构造湿空气对象
   * @param dryBulbTemperature 干球温度,单位：℃
   * @param relativeHumidity 相对湿度，单位：%；范围0~100
   * @param pressure 大气压力，单位：KPa；默认为101.325KPa
   */
  constructor(
    dryBulbTemperature: DryBulbTemperature,
    relativeHumidity: RelativeHumidity,
    pressure?: number
  )
  /**
   * ### 根据干球温度、焓值、大气压力构造湿空气对象
   * @param dryBulbTemperature 干球温度,单位：℃
   * @param enthalpy 焓值，单位：KJ/Kg
   * @param pressure 大气压力，单位：KPa；默认为101.325KPa
   */
  constructor(
    dryBulbTemperature: DryBulbTemperature,
    enthalpy: Enthalpy,
    pressure?: number
  )
  /**
   * ### 根据湿球温度、含湿量、大气压力构造湿空气对象
   * @param wetBulbTemperature 湿球温度，单位：℃
   * @param relativeHumidity 相对湿度，单位：%；范围0~100
   * @param pressure 大气压力，单位：KPa；默认为101.325KPa
   */
  constructor(
    wetBulbTemperature: WetBulbTemperature,
    relativeHumidity: RelativeHumidity,
    pressure?: number
  )
  /**
   * ### 根据湿球温度、含湿量、大气压力构造湿空气对象
   * @param wetBulbTemperature 湿球温度，单位：℃
   * @param moisture 含湿量，单位：Kg/Kg
   * @param pressure 大气压力，单位：KPa；默认为101.325KPa
   */
  constructor(
    wetBulbTemperature: WetBulbTemperature,
    moisture: Moistrue,
    pressure?: number
  )
  /**
   * ### 根据湿球温度、露点温度、大气压力构造湿空气对象
   * @param wetBulbTemperature 湿球温度，单位：℃
   * @param dewPointTemperature 露点温度，单位：℃
   * @param pressure 大气压力，单位：KPa；默认为101.325KPa
   */
  constructor(
    wetBulbTemperature: WetBulbTemperature,
    dewPointTemperature: DewPointTemperature,
    pressure?: number
  )
  /**
   * ### 根据湿球温度、含湿量、大气压力构造湿空气对象
   * @param relativeHumidity 相对湿度，单位：%；范围0~100
   * @param moistrue 含湿量，单位：Kg/Kg
   * @param pressure 大气压力，单位：KPa；默认为101.325KPa
   */
  constructor(
    relativeHumidity: RelativeHumidity,
    moistrue: Moistrue,
    pressure?: number
  )
  /**
   * ### 根据干球温度、焓值、大气压力构造湿空气对象
   * @param relativeHumidity 相对湿度，单位：%；范围0~100
   * @param enthalpy 焓值，单位：KJ/Kg
   * @param pressure 大气压力，单位：KPa；默认为101.325KPa
   */
  constructor(
    relativeHumidity: RelativeHumidity,
    enthalpy: Enthalpy,
    pressure?: number
  )
  /**
   * ### 根据含湿量、焓值、大气压力构造湿空气对象
   * @param moistrue 含湿量，单位：Kg/Kg
   * @param enthalpy 焓值，单位：KJ/Kg
   * @param pressure 大气压力，单位：KPa；默认为101.325KPa
   */
  constructor(moistrue: Moistrue, enthalpy: Enthalpy, pressure?: number)
  constructor(p1: AirState, p2: AirState, pressure: number) {
    this.pressure = pressure || Air.atm
    if (p1 instanceof DryBulbTemperature && p2 instanceof Moistrue) {
      this.t = p1.value
      this.d = p2.value
      this.h = Air.enthalpy(this.t, this.d)
      this.tw = Air.tw(this.h, this.pressure)
      this.rh = Air.rh(this.t, this.tw, this.pressure)
      this.pv = Air.pv(this.t, this.tw, this.pressure)
      this.td = Air.td(this.pv)
    } else if (
      p1 instanceof DryBulbTemperature &&
      p2 instanceof RelativeHumidity
    ) {
      this.t = p1.value
      this.rh = p2.value
      this.pv = (this.pvs * this.rh) / 100.0
      this.td = Air.td(this.pv * 1000)
      this.d = Air.moistrue(this.pv, this.pressure)
      this.h = Air.enthalpy(this.t, this.d)
      this.tw = Air.tw(this.h, this.pressure)
    } else if (
      p1 instanceof DryBulbTemperature &&
      p2 instanceof WetBulbTemperature
    ) {
      this.t = p1.value
      this.tw = p2.value
      this.pv = Air.pv(this.t, this.tw, this.pressure)
      this.rh = (this.pv / this.pvs) * 100.0
      this.td = Air.td(this.pv)
      this.d = Air.moistrue(this.pv, this.pressure)
      this.h = Air.enthalpy(this.t, this.d)
    } else if (
      p1 instanceof DryBulbTemperature &&
      p2 instanceof DewPointTemperature
    ) {
      this.t = p1.value
      this.td = p2.value
      this.pv = Air.pvs(this.td)
      this.rh = (this.pv / this.pvs) * 100.0
      this.tw = Air.tw(this.h, this.pressure)
      this.d = Air.moistrue(this.pv, this.pressure)
      this.h = Air.enthalpy(this.t, this.d)
    } else if (p1 instanceof DryBulbTemperature && p2 instanceof Enthalpy) {
      this.t = p1.value
      this.h = p2.value
      this.tw = Air.tw(this.h, this.pressure)
      this.rh = Air.rh(this.t, this.tw, this.pressure)
      this.pv = Air.pv(this.t, this.tw, this.pressure)
      this.d = Air.moistrue(this.pv, this.pressure)
      this.td = Air.td(this.pv)
    } else if (
      p1 instanceof WetBulbTemperature &&
      p2 instanceof RelativeHumidity
    ) {
      this.tw = p1.value
      this.rh = p2.value
      this.t = Air.evalTdbFromRh(this.rh, this.tw, this.pressure)
      if (this.t < this.tw) {
        throw new EvalError(
          'Wet bulb temperature is greater than dry bulb temperature'
        )
      }
      this.pv = Air.pv(this.t, this.tw, this.pressure)
      this.d = Air.moistrue(this.pv, this.pressure)
      this.h = Air.enthalpy(this.t, this.d)
      this.td = Air.td(this.pv)
    } else if (p1 instanceof WetBulbTemperature && p2 instanceof Moistrue) {
      this.tw = p1.value
      this.d = p2.value
      this.t = Air.evalTdbFromTwAndRh(this.d, this.tw, this.pressure)
      if (this.t < this.tw) {
        throw new EvalError(
          'Wet bulb temperature is greater than dry bulb temperature'
        )
      }
      this.pv = Air.pv(this.t, this.tw, this.pressure)
      this.rh = Air.rh(this.t, this.tw, this.pressure)
      this.h = Air.enthalpy(this.t, this.d)
      this.td = Air.td(this.pv)
    } else if (
      p1 instanceof WetBulbTemperature &&
      p2 instanceof DewPointTemperature
    ) {
      if (p2.value > p1.value) {
        throw new RangeError(
          'Wet bulb temperature is greater than dew point temperature'
        )
      }
      this.tw = p1.value
      this.td = p2.value
      this.d = Air.moisFromTdp(this.td)
      this.t = Air.evalTdbFromTwAndRh(this.d, this.tw, this.pressure)
      if (this.t < this.tw) {
        throw new EvalError(
          'Wet bulb temperature is greater than dry bulb temperature'
        )
      }
      this.pv = Air.pv(this.t, this.tw, this.pressure)
      this.rh = Air.rh(this.t, this.tw, this.pressure)
      this.h = Air.enthalpy(this.t, this.d)
    } else if (p1 instanceof RelativeHumidity && p2 instanceof Moistrue) {
      this.rh = p1.value
      this.d = p2.value
      this.pv = Air.pvFromMois(this.d, this.pressure)
      this.pvs = (100.0 * this.pv) / this.rh
      this.t = Air.evalTdbFromPvs(this.pvs)
      this.h = Air.enthalpy(this.t, this.d)
      this.tw = Air.tw(this.h, this.pressure)
      this.td = Air.td(this.pv)
    } else if (p1 instanceof RelativeHumidity && p2 instanceof Enthalpy) {
      this.rh = p1.value
      this.h = p2.value
      this.tw = Air.tw(this.h, this.pressure)
      this.t = Air.evalTdbFromRh(this.rh, this.tw, this.pressure)
      if (this.t < this.tw) {
        throw new EvalError(
          'Wet bulb temperature is greater than dry bulb temperature'
        )
      }
      this.pv = Air.pv(this.t, this.tw, this.pressure)
      this.d = Air.moistrue(this.pv, this.pressure)
      this.td = Air.td(this.pv)
    } else if (p1 instanceof Moistrue && p2 instanceof Enthalpy) {
      this.d = p1.value
      this.h = p2.value
      this.tw = Air.tw(this.h, this.pressure)
      this.t = Air.evalTdbFromTwAndRh(this.d, this.tw, this.pressure)
      if (this.t < this.tw) {
        throw new EvalError(
          'Wet bulb temperature is greater than dry bulb temperature'
        )
      }
      this.pv = Air.pv(this.t, this.tw, this.pressure)
      this.rh = Air.rh(this.t, this.tw, this.pressure)
      this.td = Air.td(this.pv)
    }
  }
  //#endregion

  //#region 湿空气相关的静态方法
  /**
   * ### vapor saturated pressure 饱和水蒸气分压力
   * #### 根据干球温度计算饱和水蒸气分压力；
   * @param t 干球温度，单位：℃
   * @returns 饱和水蒸气分压力，单位：kPa
   */
  static pvs(t: number) {
    // 将温度加上T0
    const ta = t + this.T0
    // 定义常数
    const c1 = -5.6745359e3,
      c2 = 6.3925247,
      c3 = -9.677843e-3,
      c4 = 6.2215701e-7,
      c5 = 2.0747825e-9,
      c6 = -9.484024e-13,
      c7 = 4.1635019,
      c8 = -5.8002206e3,
      c9 = 1.3914993,
      c10 = -4.8640239e-2,
      c11 = 4.1764768e-5,
      c12 = -1.4452093e-8,
      c13 = 6.5459673
    // 如果温度小于0.01
    if (t < 0.01) {
      // 返回计算结果
      return (
        0.001 *
        Math.exp(
          c1 / ta +
            c2 +
            c3 * ta +
            c4 * Math.pow(ta, 2) +
            c5 * Math.pow(ta, 3) +
            c6 * Math.pow(ta, 4) +
            c7 * Math.log(ta)
        )
      )
    }
    // 否则返回计算结果
    return (
      0.001 *
      Math.exp(
        c8 / ta +
          c9 +
          c10 * ta +
          c11 * Math.pow(ta, 2) +
          c12 * Math.pow(ta, 3) +
          c13 * Math.log(ta)
      )
    )
  }

  /**
   * ### vapor pressure 水蒸气分压力；
   * #### 根据干球温度和湿球温度计算水蒸气分压力；
   * @param t 干球温度，单位：℃
   * @param tw 湿球温度，单位：℃
   * @param pressure 大气压力，单位：kPa
   * @returns 水蒸气分压力，单位：kPa
   */
  static pv(t: number, tw: number, pressure = this.atm) {
    let res: number, d: number
    const pvs = this.pvs(tw)
    const ds = (pvs / (pressure - pvs)) * 0.621945
    if (tw <= 0) {
      d =
        ((2830 - 0.24 * tw) * ds - 1.006 * (t - tw)) /
        (2830 + 1.86 * t - 2.1 * tw)
      res = pressure * (d / (0.621945 + d))
    } else {
      d =
        ((2501 - 2.326 * tw) * ds - 1.006 * (t - tw)) /
        (2501 + 1.86 * t - 4.186 * tw)
      res = pressure * (d / (0.621945 + d))
    }
    return res
  }

  /**
   * 根据含湿量计算水蒸气分压力；
   * @param d 含湿量，单位：Kg/kg
   * @param pressure 大气压力，单位：kPa
   * @returns 水蒸气分压力，单位：kPa
   */
  static pvFromMois(d: number, pressure = this.atm) {
    return pressure * (d / (0.621945 + d))
  }

  /**
   * ### vapor moisture 含湿量
   * @param pvs 水蒸气分压力，单位：kPa
   * @param pressure 大气压力，单位：kPa
   * @returns 含湿量，单位：Kg/kg
   */
  static moistrue(pv: number, pressure = this.atm) {
    if (pv >= pressure) {
      throw new RangeError('pvs must be less than pressure')
    }
    return (0.621945 * pv) / (pressure - pv)
  }

  static moisfromTdbAndTw(t: number, tw: number, pressure = this.atm) {
    const pv = this.pv(t, tw, pressure)
    return this.moistrue(pv, pressure)
  }

  /**
   * ### dew point 露点温度
   * #### 根据露点温度计算含湿量；
   * @param td 露点温度，单位：℃
   * @returns 含湿量，单位：g/kg
   */
  static moisFromTdp(td: number) {
    const c1 = 9.164e-3,
      c2 = 1.446e-4,
      c3 = 1.741e-6,
      c4 = 5.195e-8
    return (
      0.001 *
      (3.703 +
        0.286 * td +
        c1 * td * td +
        c2 * Math.pow(td, 3) +
        c3 * Math.pow(td, 4) +
        c4 * Math.pow(td, 5))
    )
  }

  /**
   * ### 根据干球温度和含湿量计算湿空气焓值；
   * @param t 干球温度，单位：℃
   * @param d 含湿量，单位：Kg/kg
   * @returns 焓值，单位：kJ/kg
   */
  static enthalpy(t: number, d: number) {
    return this.Cpa * t + d * (this.latentHeat + this.Cpv * t)
  }

  /**
   * ### 根据水蒸气分压力计算露点温度；
   * @param pv 水蒸气分压力，单位：KPa
   * @returns 露点温度；单位：℃
   */
  static td(pv: number) {
    const y = Math.log(pv)
    return pv > 0.622082
      ? 6.54 +
          14.526 * y +
          0.7389 * Math.pow(y, 2) +
          0.09486 * Math.pow(y, 3) +
          0.4569 * Math.pow(pv, 0.1984)
      : 6.09 + 12.608 * y + 0.4959 * Math.pow(y, 2)
  }

  /**
   * ### 根据干球温度和湿球温度计算相对湿度；
   * @param t 干球温度，单位：℃
   * @param tw 湿球温度，单位：℃
   * @param pressure 大气压力，单位：kPa
   * @returns 相对湿度，范围：0~100
   */
  static rh(t: number, tw: number, pressure = this.atm) {
    if (tw >= t) {
      return 100
    }
    return (100 * this.pv(t, tw, pressure)) / this.pvs(t)
  }

  /**
   * ### 根据干球温度和水蒸气分压力计算湿空气密度；
   * @param t 干球温度，单位：℃
   * @param pv 水蒸气分压力，单位：KPa
   * @param pressure 大气压力，单位：kPa
   * @returns 湿空气密度，单位：Kg/m3
   */
  static density(t: number, pv: number, pressure = this.atm) {
    if (pv >= pressure) {
      throw new RangeError('pv must be less than pressure')
    }
    const ta = t + this.T0
    return (3.484 * pressure - 1.34 * pv) / ta
  }

  /**
   * ### 根据焓值计算湿球温度；
   * @param h 焓值，单位：kJ/kg
   * @param pressure 大气压力，单位：kPa
   * @returns 湿球温度；单位：℃
   */
  static tw(h: number, pressure = this.atm) {
    const acc = 5e-4
    let low = -100,
      high = 100,
      twb = NaN
    for (let i = 0; i < 1000; i++) {
      const dt = (high - low) * 0.5
      twb = low + dt
      const pvs = this.pvs(twb)
      const ds = this.moistrue(pvs, pressure)
      const hmid = this.enthalpy(twb, ds)
      if (h - hmid > 0) {
        low = twb
      } else {
        high = twb
      }
      if (Math.abs(dt) < acc) {
        break
      }
    }
    return twb
  }

  protected static evalTdbFromRh(rh: number, tw: number, pressure = this.atm) {
    let l = -100,
      h = 100,
      t = NaN
    do {
      t = (l + h) * 0.5
      const fl = this.rh(l, tw, pressure) - rh
      const ft = this.rh(t, tw, pressure) - rh
      if (fl * ft < 0) {
        h = t
      } else {
        l = t
      }
    } while (Math.abs(h - l) > Air.EPSILON)
    t = (l + h) * 0.5
    return t
  }

  protected static evalTdbFromTwAndRh(
    d: number,
    tw: number,
    pressure = this.atm
  ) {
    let l = -100,
      h = 100,
      t = NaN
    do {
      t = (l + h) * 0.5
      const fl = this.moisfromTdbAndTw(l, tw, pressure) - d
      const ft = this.moisfromTdbAndTw(t, tw, pressure) - d
      if (fl * ft < 0) {
        h = t
      } else {
        l = t
      }
    } while (Math.abs(h - l) > Air.EPSILON)
    t = (l + h) * 0.5
    return t
  }

  protected static evalTdbFromPvs(pvs: number) {
    let l = -100,
      h = 200,
      t = NaN
    do {
      t = (l + h) * 0.5
      const fl = this.pvs(l) - pvs
      const ft = this.pvs(t) - pvs
      if (fl * ft < 0) {
        h = t
      } else {
        l = t
      }
    } while (Math.abs(h - l) > Air.EPSILON)
    t = (l + h) * 0.5
    return t
  }

  /**
   * ### 根据海拔高度计算大气压力；
   * @param altitude 海拔高度，单位：m
   * @returns 大气压力，单位：kPa
   */
  static pressureFromAltitude(altitude: number) {
    // return (
    //   Air.atm *
    //   Math.pow(
    //     1 - 2.25577e-5 * altitude * (6357 / (6357 + altitude / 1000)),
    //     5.2568
    //   )
    // )
    return (
      Air.atm * (5.3788e-9 * altitude * altitude - 1.1975e-4 * altitude + 1)
    )
  }

  static densityFromAltitude(altitude: number) {
    const r0 = 1.225,
      t0 = 288.15
    if (altitude <= 11000) {
      const t = t0 - 0.0065 * altitude
      return r0 * Math.pow(t / t0, 4.25588)
    } else if (altitude <= 20000) {
      return 0.36392 * Math.exp((11000 - altitude) / 6341.62)
    } else {
      const t1 = 216.65 + 0.001 * (altitude - 20000)
      return 0.088035 * Math.pow(t1 / 216.65, -35.1632)
    }
  }

  //#endregion
}
